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ABSTRACT 

The conventional method of imposing time dependent boundary conditions for Runge- 
Kutta (RK) time advancement reduces the formal accuracy of the space-time method to first 
order locally, and second order globally, independently of the spatial operator. This counter 
intuitive result is analyzed in this paper. 

Two methods of eliminating this problem are proposed for the linear constant coefficient 
case: 1) impose the exact boundary condition only at the end of the complete RK cycle, 
2) impose consistent intermediate boundary conditions derived from the physical boundary 
condition and its derivatives. The first method, while retaining the RK accuracy in all 
cases, results in a scheme with much reduced CFL condition, rendering the RK scheme 
less attractive. The second method retains the same allowable time step as the periodic 
problem. However it is a general remedy only for the linear case. For non-linear hyperbolic 
equations the second method is effective only for for RK schemes of third order accuracy or 
less. Numerical studies are presented to verify the efficacy of each approach. 


^his research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NASI- 19480 while the authors were in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681. 


1 




Introduction 

The recent interest in long time integration is due to the need to tackle problems in areas 
such as aero-acoustics, electro-magnetics and others. This in turn necessitates working with 
higher order accurate spatial differencing operators. In many cases the time-stepping algo- 
rithm of choice is a multi-stage Runge-Kutta (RK) of temporal order of accuracy comparable 
to the spatial one. 

Several bothersome issues arise when using RK methods for long time integration. A 
principle concern is the imposition of numerical boundary conditions which retain the formal 
accuracy of the numerical method and guarantee numerical stability of the solution. For 
example (see Trefethen [1], or Carpenter et. al [2], [3]) Lax and GKS stability are not 
sufficient to assure that the there is no exponentially growing temporal error when using 
realistic grids. To alleviate this “anomaly” the use of spatial operators which have specific 
semi-discrete energy norms has been proposed [4], [5], [6], [7]. These papers have primarily 
focused on the semi-discrete form of the equations. 

The present paper, however, deals with the loss of accuracy due to the imposition of 
time dependent boundary conditions g(t), dictated by the physics of the problem. The 
conventional way of dealing with the uncertainty of what is happening at the intermediate 
stages of the RK time advancement is to impose at t + a u St, the boundary value g(t + a^St), 
where is the coefficient appropriate to the particular u th stage of the given RK algorithm. 

It will be shown that this conventional boundary condition imposition leads to a numerical 
scheme which is only first order accurate in the neighborhood of the boundary, leading to 
a global accuracy of second order only. Another approach is to treat the time-dependent 
boundary condition, g(t), as a source term in the governing partial differential equation 
(p.d.e), thereby avoiding the need to formally^ specify ijider^nediaie boundary conditions. 
However, it can be shown that procedure is equivalent to the conventional method with its 
attended problems (see Appendix A). A third natural procedure is indeed not to specify^ any 
intermediate boundary condition, but to obtain them from the inner scheme. This method 
retains the accuracy of the spatial operator, but significantly reduces the allowable time step 
for stability, rendering the scheme less attractive. 

In section 2) we analyze and pinpoint the reasons for the deterioration of the accuracy and 
provide a simple recipe’ for restoring the accuracy in the case of linear, constant coefficient 
hyperbolic system of p.d.e.’s. 

In section 3) we deal with a non-linear hyperbolic system of conservation laws. The reme- 
dies provided in Sections 2) and 3) can not be generalized to non-linear systems integrated 
by RK schemes of arbitrary order of accuracy. We show that for the RK of third order, the 
remedy of section 2) is effective. 
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2. The linear case 

To illustrate the phenomenon of loss of accuracy due to the conventional imposition of 
inflow boundary conditions, we consider the following problem: 


£ + i - 0 (» 

u(0,t) = g(t) (2) 

The semi-discretized version of equations (1) - (2) is 

dv * 

= (QV(i))i i = l,...,N;t>0 (3) 

Vo(t) = 9{t), ( 4 ) 


where V = i>, T = [u 0 > ...w/v] r is the semi-discrete approximation which converges to u(x t , t) 
at the spatial grid points Xi (for stable discretizations); and Q is the differentiation matrix 
representation of the derivative operator The specific form of Q depends on the algorithm 
used and in particular on the order of accuracy. For all finite difference operators on uniform 
grids (which suffices for the present purpose of illustration) we may write Q — j^D , where 
h is the mesh spacing. 

The demonstration of accuracy deterioration will be shown for the four-stage “classical” 
RK scheme, which is one of the most common RK time advancing schemes. For the analysis 
to make sense we assume that the differentiation matrix Q, is at least of fourth-order accuracy 
up to the boundary. It should be noted, however, that this illustration could be carried out 
for any RK algorithm. 

The above mentioned four-stage integration, together with the conventionally imposed 
boundary conditions, is implemented as follows: 


v! 

= 

+ 

kw), 

1 < i < N 

(5) 

v o 

= g{t 

+ 

*, 
2 ’ 



vf 

= V? 

+ 

Vk 

VI 
* «** 

VI 

i-H 

(6) 

v 2 0 

= g(t 

+ 

2 ’ 



Vi 

= V? 

+ 

A {DV% 

1 < i < N 

(7) 
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V 0 3 = g(t + St) 

u? +1 = v? + ^[{DV n ) i + 2(DV 1 ) i + 2(DV 2 ) i + {DV 3 ) i } 1 < t < N (8) 

0 

u o +1 = s(* + tit) 

where A = 

Equations (5) - (8) take the semi-discrete variable v,-(f), from the time level n, to u,(t + £t) 
at time level n + 1 . 

For the purpose of analysis, the above system is rewritten in the following form, again 
with V = [uoi ••• u w] • 


V 1 = K n + £ DV n + G°e 0 (9) 

Li 

V 2 = K n + ^DV l + G 1 e 0 (10) 

V z = v n + \DV 2 + G 2 e 0 (11) 

yn+l _ yn + K D yO + W yl + 2D y2 + D y3 J + q3 6q ( 12 ) 

6 

where eo = [1,0, ...0] , and 

G° = g(t + ~)-vS- \{DV") „ (13) 

G' = g(t + |) - v 0 "-l(W)o (14) 

G 2 = 9 (i + St) - «J - \(DV 2 )o (15) 

G 3 = j(< + St) - vj - 1[W“ + 2W + 2DV* + Z?V 3 ]o (16) 

0 


Substitution of (9) into (10) and the result into (11), etc. leads to the following expression 
for l/ n+1 : 


l/” +1 = [/ + ^ D 2 + ^D 3 + £ D'\V « 

+ g°[^d + £» 2 + 

+ G '[? D + T D ^ e ° 

O 0 

+ G 2 [jD)e° 

+ G 3 e° (17) 
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To determine the formal accuracy of (17) we substitute the exact solution u(xi,t), (pro- 
jected at the points Xi), for v”. Under the assumption on the order of the differentiation 
matrix Q, it is clear that D k V n in equation (17) can be replaced, up to fourth-order accuracy, 
by h k [(-\) k - g at the points (x,-,t). Equation (17) then becomes 


1/"+’ = u (t + St) + 0([£<] 5 , [&r] 5 ) 

+ <?{\d + X jD> + ^] e ° 

+ G'&D + 

O 0 

+ G 3 (jD]e° 

+ G 3 e° (18) 

Applying the same procedures to equations (13) - (16) we get: 


G° — g(t + — ) - g(t) - yfl('(<) 

G' = S(t + | )-»(<)- { -fm - - £<&, G» 

G ! = g(t + St) - g(i) - (St)g'(t) - ^ g"(t ) - ^g'"(t) - A<, G' - ^4, G° 

G 3 = g(t + St) - g(t) - (St)g'(t) - fg"(t ) - ^V'(i) - ®V"«) 

jl . (^) 2 j 2 | (^) 3 j 3 l /~(0 

l^doo + “g” d oo + 

- [g^oo + ^~^oo] 

- jdlo G 3 (19) 

where d^, = (De 0 ) 0 , = (D 2 eo)o, do 0 = (D 3 e o)o, and (D k e o)o is the first element of 

the vector D k e o; 1 < fc < 3. A Taylor series expansion of equation (19) clearly shows that 
G° — 0([df] 2 ); as it does for G 1 , G 2 and G 3 , for arbitrary A. In addition, it can be shown 
that the vectors D k e o (1 < k < 3) are linearly independent. Substituting equation (19) into 
equation (18) yields the expression 


V ntl - u(t + St) = [St\‘^-(0 ! e o )g" + 0([«] 3 ) 

yt> 


( 20 ) 
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Rearranging equation (20) we have 


— > = Sx { -^(D\ a )g" + 0([«i] 2 ) = 0(6x) (21) 

which means that the RK approximation is locally (near the boundary) only first order 
accurate. The locality is assured by the finite support of D for explicit schemes, and by geo- 
metrically decaying coefficients in space for compact schemes. However, the overall accuracy 
in the L; n f norm cannot exceed second order (see Gustafsson [8]). 

The remedy for the dilemma posed by equation (21) suggests itself when one examines 
equations (17) and (18). We note that if we set G° = G 1 = G 2 = 0, then V n+l = 
u(t + St) + 0([<5f] 5 ) + e 0 G 3 . But, with G° = G l = G 2 = 0, then G 3 = 0([6f] 5 ), and we 
have the correct order for V rn+1 . To achieve these identities we specifically use in equations 
(5)-(8) the following expressions for the intermediate boundary conditions: 


= s(0 + js'W 

(22) 

= 9(0 + + ^9"(0 

(23) 

— g(t ) + (8t)g (t) + ^ 9 (t) + ^ 9 (0 

(24) 


Note that equations (22 - 24) are precisely the intermediate boundary conditions obtained by 
numerically integrating ^ = g"'(t ) with the classical fourth order RK scheme. Specifically, 
replacing the third order equation in u with the system of first order equations for the 
unknown functions Vo, vt 0 and vtto, 


d(V) o 

dt 

d(vt) 0 

dt 

d(vtt ) o 

dt 


(vt) o 
(utt)o 
<?"'(*) 


(25) 

(26) 
(27) 


on the boundary, and integrating with the standard fourth order RK scheme, yields precisely 
the intermediate boundary conditions proposed in equation (24). Solving the third order 
o.d.e (27) on the boundary, is a remedy that can be applied to any four stage fourth order 
RK scheme, and thus provides a simple and general means of implementing the correct inter- 
mediate boundary conditions. For three stage third order RK schemes, solving a It 2 — 9 (0 
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on the boundary provides a general technique for imposing the correct boundary conditions 
on the linear problem. 

At this point we comment on why the above predicted phenomena has not been observed 
previously by the practitioners of higher order methods. From equation (21) we see that 
the leading error coefficient in the expression is proportional to (A) 4 £x. For example, if D 
represents a fourth order central difference operator with fourth order boundary closures, 
and it = then the error at the point next to the boundary becomes 


*(184320 - 1520A)A 4 6x 2 (+172800 + 23040A - 190A 2 )A 4 6x 3 

Ut + Ux ~ 165888 + 79626240 + 79626240 + • • -( 2 ) 

Thus if A << 1 for a given 6x it is possible that \ 4 Sx < (£f) 4 . However, if one refines the 
grid, or conversely, if one runs the computation at the allowable stable CFL (A = 0(1)), 
then the first order error becomes apparent. A detailed study of this effect is now presented. 

We begin by showing the the results from a grid refinement study using three high-order 
central difference schemes; 1) (3, 3-4-3, 3): a fourth order explicit interior with two third-order 
explicit boundary points at each end of the domain, 2) (5 2 ,5 2 — 6 — 5 2 ,5 2 )[2]: a sixth-order 
compact interior, with two fifth-order boundary closures at each end of the domain, and 3) 

(5, 5, 5, 5, 5, 5-6-5, 5, 5, 5, 5, 5)[6]: a sixth-order explicit, with six fifth-order boundary points at 
each end of the domain. (See Carpenter et al[2] for detail of the high-order nomenclature). 

All schemes are GKS and time stable for the scalar hyperbolic equation. In addition, when 
used in conjunction with the SAT [3] boundary procedure, the semi-discrete form of scheme 
3) can be shown to be time stable for the constant coefficient hyperbolic system. An explicit, 
a compact implicit, and a scheme which satisfies the summation by parts energy norm were 
chosen to illustrate the generality of the problem, as well as the remedy. 

The model problem used to test the schemes is the scalar hyperbolic equation 


u u t/u 

dt ^ dx 

= 0, 0 < x < 1, t > 0 

(29) 

u(0,t) 

= sin 2n(— t), t > 0 

(30) 

u(x,0) = 

sin 27r(x), 0 < x < 1, 

(31) 


The exact solution is 
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u(x,t) = sin27T {x-t), 0 < x < 1, <>0 (32) 

In all cases the solution was advanced in time using the classical fourth-order RK scheme 
with boundary conditions imposed as in equations (5) - (8). 

Upon grid refinement with a vanishingly small CFL, all schemes recover their theoretical 
accuracy of fourth-order, sixth-order and sixth-order, respectively. Table (I. a) shows that 
all schemes converge at a fourth-order rate for moderate resolutions, but gradually degrade 
to an asymptotic rate of 2.5. This trend is representative of high-order explicit or compact 
spatial operators advanced in time with a fourth-order RK scheme, with or without the SAT 
boundary procedure. 



(3, 3-4-3, 3) 


(5 2 ,5 2 - 6 - 5 2 ,5 2 ) 


(5, 5, 5, 5, 5-6-) 



CFL = 2.0 

Conv 

CFL = 1.4 

Conv 

CFL = 1.5 

Conv 

Grid 

log L 2 

Rate 

log L 2 

Rate 

log L 2 

Rate 

41 

-2.371 


-3.033 


-1.537 


81 

-3.570 

3.98 

-4.242 

4.02 

-2.221 

2.27 

161 

-4.678 

3.68 

-5.450 

4.02 

-2.956 

2.44 

321 

-5.593 

3.04 

-6.635 

3.94 

-3.704 

2.48 

641 

-6.380 

2.61 

-7.711 

3.57 

-4.455 

2.49 

1281 

-7.138 

2.52 

-8.586 

2.91 

-5.208 

2.50 


Table I.a: Conventional imposition of boundary condition u(0, t) — g{t ), for three 

high-order finite difference schemes. 


Two possible remedies have been suggested to rectify the loss of accuracy for the linear 
problem: 1) do not impose intermediate physical boundary conditions, and 2) impose con- 
sistent intermediate boundary conditions derived from the physical boundary conditions and 
their derivatives. (Or alternatively, solve the related system of equations on the boundary). 

Not imposing intermediate physical boundary conditions results in a fully discrete scheme 
which is formally fourth order accurate. (By definition a fourth order scheme in space and 
time will remain fourth order in the absence of boundary conditions). A problem with this 
remedy, however, is that the stability of the scheme is greatly reduced. When using the RK4 
scheme with various finite difference operators, at least of a factor of two (and often much 
more) decrease in CFL was observed. 

The alternative remedy is to impose consistent intermediate boundary conditions, derived 
from the physical boundary conditions and its derivatives. For the scalar advection defined 
by equations (29) - (30) it is sufficient to solve the derivative boundary conditions described 
by the system of equations (27). 
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Table (I.b) shows the results of a grid refinement study using the derivative form of the 
boundary condition. 



(3, 3-4-3, 3) 


(5 2 ,5 2 — 6 — 5 2 ,5 2 ) 


(5, 5, 5, 5, 5-6-) 



CFL = 2.0 

Conv 

CFL = 1.4 

Conv 

CFL = 1.5 

Conv 

Grid 

log L 2 

Rate 

log L 2 

Rate 

log L 2 

Rate 

41 

-2.394 


-3.090 


-2.490 


81 

-3.613 

4.05 

-4.282 

3.97 

-3.767 

4.24 

161 

-4.817 

4.00 

-5.486 

3.99 

-5.076 

4.35 

321 

-6.019 

3.99 

-6.687 

3.99 

-6.377 

4.32 

641 

-7.222 

3.99 

-7.891 

4.00 

-7.655 

4.25 

1281 

-8.425 

4.00 

-9.099 

4.01 

-8.911 

4.17 


Table I.b: New physical boundary condition imposed as - = g " (<), for three 

high-order finite difference schemes. 

Fourth-order temporal accuracy is recovered for all methods with this approach. In addition, 
the same maximum CFL was achieved in all cases, as was possible with the conventional 
boundary conditions. Similar results have been obtained for three stage third order and six 
stage fifth order RK schemes as well. 

Section 3. Non-Linear Hyperbolic System 

We have shown that the traditional imposition of time-dependent boundary conditions 
causes a degradation of the formal accuracy to first-order. We have also shown how to elim- 
inate the problem, given the exact boundary condition and its derivatives on the boundary. 
We shall now show that boundary treatment similar to that resulting from the linear analysis 
of section 2, is also valid for third order RK schemes, even for the case of a system of non- 
linear conservation laws. The procedure does not obviously generalize to higher order RK 
integration, however. (The only technique available that does directly extend to non-linear 
hyperbolic systems is not imposing intermediate physical boundary conditions.) 

Consider the system 


f = ^O) 9 £=r x m,0< x <l;t>0 (33) 

0(0, t) = g(t) (34) 

where A(U) is the Jacobian of f(U) with respect to U . We perform the integration using the 
third order RK scheme attributed to Heun. The boundary conditions at the intermediate 
time levels are obtained from solving the boundary equation = g"(t) with Heun’s 

method. 
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= ®r + j(D/(3")i, 

1 < i < N 

(35) 

^0 

+ 

•Vi 

II 



fi? 

= ST + y P/V)]i 

1 < i < N 

(36) 

fig 

- s(0 + yiflO + 



u" +1 

= W + 7^; + jA[D/V)li 1 < * < 

4 4 4 

(37) 

i7 n+1 

u 0 

= g(t + $t) 




Following the notation of equations (9) - (12) we have* 


$ = n? + ^[Df(^)]i + G°e o (38) 

vf = v? + -^-[Dfiv 1 )] i + G 1 e 0 (39) 

+ 1$ + ^[Df{^)]i + & eo (40) 

where, it is clear from equation (35) - (37) that here: 

<5° = m + js'w - <5 - jP/tojo ( 4 >) 

s = m + - ywV)]° (- 12 ) 

S 1 = s(j + ft) - - yj - jA[D/(t?)l„ (43) 


As in the linear case described by equation (18), G 2 will be given by a linear combination 
of G° and G 1 , plus terms proportional to (6f) 4 . We now show from equations (41) and (42) 
that G° = 0, and G 1 = 0([<5t] 3 ), and thus the overall accuracy of G 2 is 0([<5t] 3 ). 

As in section 2, in order to obtain the truncation error we substitute u(x i,t) for uf . The 
vector G° immediately follows as 

G°(«) = g{t) + -jj"if(*) — (u(t) + — ,/r(«) )o (44) 

= 9(0 + fm - (3(0 + j S.)o = 0 (45) 
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and in a similar vain, 


<?' = m + - ( am - ^£i/( «(() + |/;(s)]}„ (4 6) 

= if(0 + + g ^ ?'(<) - { “ - ^-^[/W)) + + 0([«] 2 )]}o 

(47) 

but f s u t = v4(u)u ( = f t and £[f t ] = ^(/*) = u tt . Thus 


G 1 

G 1 


m + + 

0 ([^] 3 ) 


2(sty 


m - { $ 


28t^ 

-T“‘ + 


2(6tf- 


i^tt}o + G([<5t] 3 ) 


(48) 


The consequence of equations (45) and (48), substituted into equation (40) with t/" — > u(x, t) 
is that the error is proportional to 


\vr^ - uijt + st ) i 

St 


0([St] 2 y, i = 0,1 


m 


(49) 


near the boundary, (m finite), and proportional to 0([6£] 3 ) in the interior. In other words 
, the boundary conditions prescribed in equations (35) - (37) give us the required overall 
accuracy. Unfortunately, this procedure does not seem to generalize to RK integrations of 
order higher than three for the case of non-linear systems. The main reason is that beyond 
u t t, we will need to use the “Jacobian of the Jacobian” which can not be related simply to 
the temporal derivatives of the vector u(x,t). 


4. Conclusion 

We have shown that the imposition of inflow boundary conditions at the intermediate 
steps of Runge-Kutta algorithms for hyperbolic P.D.E.’s has to be done in a counter-intuitive 
manner, if one is to preserve the overall accuracy of the scheme. The conventional (or 
"natural”) way of assigning at level n + ati, (0 < a, < 1) the value g(t + a,St) degrades 
the scheme to be of first order accuracy near the boundary and of second order accuracy 
overall. We presented ways to deal with the linear case of general order and with systems of 
conservation laws in the case of third order RK integration. Much work remains to be done 
for the non-linear case and linear problems with variable coefficients. 
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Appendix A 

Once the spatial operator has been chosen, a P.D.E. becomes a system of O.D.E. s, plus 
a boundary term. If the boundary term is treated as a source term, then entire system can 
then be treated by conventional techniques. We show that this technique suffers from the 
same loss of accuracy at the boundary as was discussed earlier. 

We start with the governing equations 


du du 
dt dx 


u(0, t ) 


0 0 < x < 1 , t > 0 

9{t) 


The semi-discretized version of (AI. 1) - (AI. 2) is 


(AI. 1) 
(AI. 2) 


•ti = i (Dv), > = 1 JV; * > 0 (AI- 3) 

at h 

Vo (t) = g(t) (A'- 4 ) 


where V = u, T ,i = 0, ..., N. Define the reduced identity matrix P by 


P = 



1 


where the matrix P is dimension (N+1,N-|-1). Similarly, eo = [1,0,...,0] , and has dimension 
(N+l). Next define the new variables V = PV, D = PDP and qo — PDe o. Now replace 
equation (AI. 1) by the fully discrete source equations: 


V? 

= V? 


+ S«(«o)il 1 < i < N 

(AI. 5) 

V? 

= V? 

+ 5 l(bv 

+ + y)(9o)j] 

(AI. 6) 

V? 

II 

+ M(W)’ 

+ 9{t + _ 2“)(4'°) j] 

(AI. 7) 

v n+1 

— V n 

r J 

+ g[(BV n ). 

i + 2(DV l )j + 2 (DV 2 ), + (DV 3 ),} 



+ j[<7(*) +4^(f + — ) -1- g(t + ^)](9o)i 

0 X, 

(AI. 8) 
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Rewriting equation (AI. 5) in terms of V and D one obtains: 


(PV'), = (PV n ), + j(! (PDPV’ '), + g(t)PD(e „),] 1 < j < N 

= (PV"), + \(PDV’% = />[(^“) + ±(DV")), (AI. 9) 

It is clear that the reduced vector PV 1 is identical to the vector at time level 1 obtained 
from the conventional imposition of boundary conditions for i <j < N [see equation (5)]. 
Noting this equivalence, equation (AI. 6) can now be interpreted as 


P(V 2 h = P(V"),A \{(PDPV'), + 9 (t + j)PD(e 0 ),} l<j<N 

= P(V")i + \(PDV'), = P[(V n ) + j(OV')],] (AI. 10) 

The reduced vector PV 2 is identical to that obtained with the conventional boundary condi- 
tions for 1 < j < N. This procedure can be used to show that each stage of the conventional 
boundary procedure is equivalent to that obtained from solving equations (AI. 5) - (AI. 8). 
Thus, the entire procedures are equivalent. 
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